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Abstract 



In this paper wc propose and analyse composite Filon-Clcnshaw-Curtis quadrature rules 
for integrals of the form ij. (/,<?) := f(x) exp(ikg(x))dx, where k > 0, / may have 
integrable singularities and g may have stationary points. Our composite rule is defined on a 
mesh with M subintervals and requires MN+1 evaluations of /. It satisfies an error estimate 
of the form C^k~ r M~ N ~ 1+r , where r is determined by the strength of any singularity in / 
^ c*j and the order of any stationary points in g and Cn is a constant which is independent of k 

and M, but depends on N. The regularity requirements on / and g are explicit in the error 
estimates. For fixed k, the rate of convergence of the rule as M — > oo is the same as would be 
obtained if / was smooth. Moreover, the quadrature error decays at least as fast as k — > oo 
as does the original integral it"' 6 '' (/, g). For the case of nonlinear oscillators g, the algorithm 
requires the evaluation of <7 _1 at non-stationary points. Numerical results demonstrate the 
sharpness of the theory. An application to the implementation of boundary integral methods 
qq \ for the high-frequency Hclmholtz equation is given. 

<N 

£NJ ' Keywords Oscillatory integrals, Clenshaw-Curtis quadrature, Integrable singularities, Station- 

£ — ■ ary points, Graded meshes. 
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1 Introduction 

Oscillatory integrals of the form 

I [ k ,b] (f>9) ■= I" f(x)exp(ikg(x)) dx (1.1) 

J a 



where / € L [a, b] and k > regularly appear in applications. If / and g are smooth and g' 
does not vanish then l]^' b \f,g) decays with at least 0(k~ 1 ) as k — > oo. The decay is faster than 
0(k~ v ) if / and some of its derivatives vanish at both end-points a, b, but is generally slower if 
/ has a singularity or if g has a stationary point in [a, b]. In practice one may be interested in 
computing (jl.ip efficiently and to controllable accuracy for a range of values of k and for quite 
general / and g. The purpose of this paper is to provide stable quadrature rules for this task and 
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to prove error estimates demonstrating the quality of the methods. We also demonstrate the 
efficiency of our rules by applying them to an example coming from boundary integral methods 
for high-frequency wave scattering. 

Restricting first to the case g(x) = x, a recent paper [5] studied the convergence of Filon- 
Clenshaw-Curtis (FCC) rules for computing the special case of ([Lip : 

4 If) := j /(x)exp(i*x) dx . (1.2) 

These rules (denoted by Ik,N(f)), approximate (jl.2p (when k > 1/2) by replacing / by its poly- 
nomial interpolant of degree N at the Clenshaw-Curtis (or Chebyshev) points t.-jv = cos(jN/ir), 
j = 0,...,N. (When k < 1/2, standard Clenshaw Curtis rules are used instead.) A stability 
theory is given in [5] and a slight extension of the error estimates given in [5] (see Theorem 12.11 
below) shows that, for r G [0,2] and m > max{i, p(r)}, 

\h(f) - hAf)\ < C[t) \\M\ H ™, N>1, (1.3) 

where f c (8) = f(cos0), p(r) = r, r G [0, 1], p[r) = 5r/2 - 3/2, r G [1, 2] and || • \\ H m denotes the 
norm of the Sobolev space of order m on [— vr, tt]. Thus fast convergence of the rule with respect 
to N and decay of the error with order up to 0(k~ 2 ) is obtained if / is sufficiently regular. 

However the convergence rate of the FCC rule is significantly impaired when / has one or 
more (integrable) singularities. Thus in this paper we consider composite rules for (jl.ip (first for 
g(x) = x), obtained by subdividing [a, b] into a mesh with M subintervals, chosen so that any 
singular points of / coincide with mesh points. We then construct a composite rule which uses 
the FCC rule on each mesh subinterval not containing the singularities, and, on subintervals 
containing the singularities, either zero or a very simple two-point rule is used, depending on 
the strength of the singularity. (See Section [3.1 1 for precise description of the algorithm.) To 
give a flavour of our results, we show, for example, that if / has a singularity of form \x — xq\^ 
for xq G [a, b] and (3 G (—1, 1)\{0}, then with suitable mesh refinement near xo, our rules have 
error E(f) which satisfies the estimate (see Theorem 13. 6p : 

where r G [0, 1 + /3] and the norm on / is an appropriate weighted norm which takes into account 
the singularity at xq . The estimate (|1.4p decays at least as fast with k as does the corresponding 
integral (jl.l|> . since the latter decays in general with 0(k~ r ) where r = min{l + /3, 1} - see 
Lemmas 13.11 - 13.31 

In order to prove (ll.4p (and its generalisations), in this paper a non-trivial extension of the 
estimate (|1.3p (quoted from [5]) is first obtained in §|2j Since the error estimate (|1.3p depends on 
the regularity of / through the norm of f c this estimate does not provide the correct scaling with 
respect to h when it is transported to an interval of size h. Therefore in Theorem 12. 51 we prove a 
variant of (jl.3p . where ||/ c ||_f/ m is replaced by the Chebyshev weighted norm of/( m ). This new 
estimate has the correct scaling behaviour, as is shown in Theorem 12.61 In §|3] we obtain the 
error analysis for the composite FCC rule applied to (|l.ip with g{x) = x when / has integrable 
singularities at a finite set of points, in particular obtaining error estimates of the form (ll.4p . In 
2] we further extend to the case where g may have a finite number of stationary points in [a, b]. 
The latter case can be reduced to that studied in §[3] provided we assume that the inverse of g 
is known (or is evaluated numerically) on subintervals between stationary points, and indeed 
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the action of g^ 1 is required for the implementation of the algorithm. In £j5]we give numerical 
experiments, utilising the public domain code [3j which indicate that as M increases, the error 
decays with 0{M~ N ~ l ), provided the parameters of the mesh are appropriately chosen, relative 
to the regularity of / and stationary points in g. Moreover, when k increases the error decays 
roughly with (D(k~ r ), where r £ [0, 1 + 0\ indicating the sharpness of our theory. The numerical 
experiments also indicate that applying the composite FCC on a graded mesh to integrals with 
singularities yields much more accurate results than applying FCC globally, and using the same 
number of integrand evaluations. 

In this paper we restrict our error estimates to the case of k > for convenience only; the 
rules also work well for all k £ M. and the error estimates can be easily extended to that case 
(see, e.g., [5j Corollary 2.3]). As is also shown in [5], the FCC rules for (|1.2p have a stable 
implementation for all k and N which, via FFT, costs 0(N log N) operations. The composite 
rules presented here require the evaluation of / at MN + 1 points. 

Although oscillatory integration is well-studied in the classical literature, some problems 
of interest to numerical analysts (even in ID) still remain unsolved today. Thus this field 
has enjoyed a recent upsurge of interest, partly because of its importance in wave scattering 
applications. (See [2] and [5] for some more detailed historical remarks.) In particular, the 
construction and analysis of Filon-type methods has been examined in Iserles |10l lllj . Iserles 
and N0rset [12], Olver [16], Xiang [21] and Huybrechs and Olver [9]. (Other related methods 
include those of Levin- type [HI [T71 [18] and those using numerical steepest descent [Sj.) In all 
these references, however, the analysis concentrates on accelerating the convergence as k — > oo, 
generally assuming either that / is sufficiently regular or / has a particular type of singularity 
so that the moments (i.e. integrals (jl.lj) where / is replaced by polynomials) can be computed 
using special functions. By contrast, we propose Filon-type method for computing (jl.ip where 
/ has algebraic singularities and g may have stationary points and where the moments can 
be obtained readily. Our method converges super algebraically with respect to the number of 
quadrature points for any strength of singularity, provided the parameters of the mesh are chosen 
appropriately, and also converges with respect to k at least as fast the integral itself converges 
to zero as k — > oo. Our error estimates explicitly indicate the regularity requirements on / and 
g. Other papers [15] and [5] provide analogous estimates for pure (non-composite) Filon rules, 
where / and g are sufficiently regular (and g' ^ ). But apart from these we know of no other 
contributions in this direction. 

Finally we mention that our methods add something to traditional asymptotic methods. 
The method of stationary phase produces an accurate approximation to an integral if k is 
sufficiently large, whereas our methods work for all k and are superalgebraically convergent 
with respect to the number of function evaluations. Our methods also yield a relative error 
which is superalgeraically convergent uniformly in k and may indeed even decay with k. 

2 The Basic Filon Clenshaw-Curtis Rule 

In this and the next section we will consider only the linear oscillator g(x) = x in (jl.ip . (See £JH 
for the case of nonlinear g.) Also, we introduce the notation 




(2.1) 



When [a,b] = [—1, 1] we will denote (|2.ip simply as Ik(f)- 
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2.1 Integrals over the fundamental interval [—1,1] 

The FCC rule in its simplest form approximates Ik(f), by replacing / by its algebraic polynomial 
interpolant QnI at the Clenshaw-Curtis points tj t N '■= cos(jn/N) , j = 0, . . . , N where N > 1. 
Then for k > 1/2, the rule is 

1 N 

hMf)--= / (Q N f)(x)exp(ikx) dx = ^2" a n>N (f)"n(k) , (2.2) 

n=0 

where, for n > 0, oj n (k) := ^ T n (x) exp(i/cx) dx , T n (x) = cos(n arccos(x)) is the nth Cheby- 
shev polynomial, and 

2 ^ 'n7r 

a n , N (f) = JrY," cos Vjf) f { tj ' N ) ' ™ = 0,..-,iV. (2.3) 

j=0 

The notation means that the first and last terms in the sum are multiplied by 1/2. 

When < k < 1/2 the integrand in (|2.ip is non-oscillatory and we then apply the standard 
Clenshaw-Curtis rule: 



1 N 

"n,Ar(/fe)w n (0) , fk{x) ■= f(x) exp(ifcx) . (2.4) 

1 „_ n 



AT 

hMf) ■-- 

n=0 

We point out that w n (0) = 2/(1 — n 2 ) if n is even, otherwise. For k > 1/2, the computation of 
oo n (k) turns out to be more delicate. However in [5] a stable and efficient scheme for computing 
these weights is presented. After an initial application of the discrete cosine transform (via 
FFT, costing 0(N\ogN) operations), the rule (|2.2p - (|2.4p can then be applied to any / in an 
additional O(N) operations - see [5] for more detail. In the error analysis in this section, we 
shall make use of the Sobolev space H m of 2ir— periodic functions with the norm 



M 2 Hm := \m\ 2 + E H 2 "V(«)| 2 , £M : = T- r V^)exp(-iu#) dO. (2.5) 
For 4> G H° = L 2 [— 7r, n] and any J > 0, we introduce the truncated Fourier series: 



J 

(Sj<l>)(e)= Y, WexpM) 
fl=—J 

which converges to <f> G H° as J — > oo. In fact when <f> is even, we have 

^ ^ /*7T 

(5j0) (0) := ^(0) + 2 V ^(n) cos(u6») , where ka) = - (f)(0) cosU6)d6 , (2.6) 

where the second sum is void if J = 0. 

Introducing the notation f c (0) = /(cos 6), and p{r) = r, r G [0, 1], and p(r) = 5r/2 — 3/2, 
r G [1, 2], the following theorem is then a minor extension of Theorem 2.2]. 

Theorem 2.1 There exists a constant C > suc/t i/tai, /or a// r G [0,2] and aZZ integers 
m > max{l/2, p(r)}, i/ie estimate 

( 1 \ r V 1 \ m ~ p(r) 

|4(/)-4,at(/)| < C(-J l-j ll/dlflm, fc>l, iV>l (2.7) 

ZioWs u>/ien / c G i7 m . 
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Proof. In [5] the estimate (\2.7\i was obtained for k > 1/2, for r = 0, 1,2. Hence for any 
9 € [0, 1] we can write 

Mf) - hMf)\ = Mf) - hAf)\ 6 Mf) - hAft' 6 ■ (2-8) 

Then using estimate (j2.7|) with r = 1 (respectively r = 2) to bound the first (respectively second) 
factor on the right hand-side of (12, 8p we obtain 

\ 2-9 / ]_ \ m-7/2+56»/2 



|Ife(/)-4,iv(/)| < C 1-1 1-1 \\f c \\ H m . 

Then setting 9 = 2 — r we obtain (|2.T[) for r 6 [1,2] and for > 1. An even simpler interpolation 
argument obtains the estimate for r € [0, 1]. For k < 1/2, i.e., for the classical Clenshaw-Curtis 
rule, the case r = follows by the same arguments used in [5j Theorem 2.2] to prove (|2.7p . For 
r E (0, 2) the result is obvious since k is bounded. □ 

Theorem 12.11 ensures arbitrarily high convergence for the FCC rule as N — > oo, provided / 
is sufficiently smooth. When / is not smooth it is better to apply the FCC rule in a composite 
fashion on meshes graded suitably towards the singular point(s). These composite rules then 
typically have fixed N and converge as the subinterval size shrinks to zero. In order to obtain 
good error estimates for the composite rules we need to modify the error estimate in Theorem 
12.11 so that derivatives of / rather than derivatives of f c appear in the bound. This will be done 
in Theorem 12. 51 which in turn is used to obtain Theorem 12.61 showing how the error of the FCC 
rule, when applied on an arbitrary interval, depends on the length of the interval. In order to 
prove Theorem 12.51 we first need two lemmas. 

Lemma 2.2 Let f be such that (/') c 6 L l [—ir,7r]. Then 

/c(A*) = ^f(70c(M-l)-(7)c(M + l)l , for MT^O. (2.9) 
Proof. Since f c is even we use (|2.6p and integrate by parts to obtain 

7c(M) = - [ W f c (9)cos(fi9)d9 = -— [ T (f c )'(9) S m(fi9)d9 

7T Jo ""A* JO 

= — [ f'(cos9)sm9sm(fi9)d9 
^ Jo 

= -J— (/') c W [cos(0* - m - C ™((V + d9 , 
and the result follows. □ 



Using Lemma 12.21 we now estimate the error in the truncated Fourier cosine series of f c . 
Lemma 2.3 For all < m < N + 1 there exist constants cr m ^ > such that 

\\(l-S N )f c \\ H m < 0-m,N\\(f im) )c\\H° ■ (2-10) 

Proof. Since Sn is the orthogonal projection of H° onto span{exp(ij^) : < \j\ < N}, the 
result is trivial for m = 0. So let us assume now that m > 1. Since f c is even, from (|2.6p we 
have, for all J > 0, 

(I-Sj)f c = 2 fc^)cos(fi9), and - S j) f c f Hm = 2 ^ » 2m \%(»)f \ 

n>j+i fi>J+i 
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Then, using Lemma 12.21 we obtain 



{i-sj)f c f Hm < \ e ^ m - 2 



(/')c(m - 1) - (/OcCm + 1) 



< E ^ 2m ~ 2 (/Oc(m-i) + E ^ 2m ~ 2 (/OcOi+i 
= E(m+ i) 2m " 2 |(rTc(M)| 2 + E (^-i) 2m_2 |(7%(M) 

/Lt>J fl>J+2 

< 2E^+i) 2m - 2 |(K(^)| 2 • 



(2.11) 



Hence, in the case J > 1, we have, 



\\(I-Sj)f c \\ 2 Hm < 2 



J + l 
J 

J + l 
J 



2m-2 



E« 

li>J 



2m-2 



2m-2 



(2.12) 



Using this identity to — 1 times (recalling that m < iV + 1), we obtain 
||(/-5jv)/c||fl™ < 

||(/-5 iV - m+1 )(/( m - 1 )) c || Jfl 



JV + l\ m -V iV \ m " 2 /iV-m + 3 



which we write as 



N - 1 



iV - m + 2 



(J — S]si)fc\\H m < CTm,iV || (J — <5jV-m+l)(/ 



m— 1)^ 



(2.13) 



Now, if m < iV+1, we can use (I2.12h one more time and use the fact that S^-m is an orthogonal 
projection on H° to obtain the required result. On the other hand if m = N + 1 we can use the 
fact that 

IKz-SoX/^y^ < ||(/ (m) ) c ||tfo , 

which is easily obtained from (I2.1ip . to deduce (I2.10|) . □ 



Remark 2.4. To estimate the constants (T m ,N> note that in each pair of terms in the product, 
we can cancel the denominator in the left-hand term with the numerator in the right-hand term 
to obtain, for m > 1, 



,771—1 



N(N — 1)...(N -m + 2) 



m—l 

n 



N + l 
N + l-jJ ' 



(with the product being interpreted as 1 when m = 1). Thus for fixed m > 1, a m> N -+ 1 as 
N — > oo. Moreover letting m grown with N (for example m = N + 1), we have 



(AT + 1)^ JV N / 1 

<TJV+1 ' JV = m = ^vTl 1 + iv 



N e N+l 



where the last relation is obtained using Stirling's formula. Since in this paper we will use fixed 
order methods (i.e. N fixed) and obtain convergence for composite methods as the mesh size 
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shrinks, the growth of cr/y+i jv is n °t of essential importance to us here. However if we wanted 
to use hp quadrature then this growth would be important and would need to be cancelled by 
suitable decay of the derivatives of / in order to obtain convergence. Estimates of this type are 

in mn. □ 



Now in Theorem 12.5 1 below we will obtain the analogue of Theorem 12 .11 but with (appropriate 
weighted norms of) derivatives of /, rather than f c on the right-hand side. For any integer m > 0, 
and a function / defined on [a, b], we introduce the weighted seminorm 



l/<™>MI 2 1 1/2 



ma *» := 17. v(t-«)V-H ' <2 ' 14) 

When [a, b] = [—1, 1] we just write | • \jjm and we note that 

= {jT|(/ (ro) )e(*)| a d0 } 1/2 = V^||(/ (m) ) c Lo- (2.15) 

Theorem 2.5 Let r G [0, 2] and < m < N + 1. There exist constants a' m N such that 

|4(/)-4,v(/)| < <Jv(t) f^J l/V™, (2-16) 

when \/\h™ < oo. Moreover a' m N = Co~ mt N with C independent ofm,N. 

Proof. Note that if \f\n™ < oo then / G L 2 [— 1,1] and so we can define an algebraic 
polynomial p of degree N by 

v 

p(x) = / c (0) + 2 ^ f c (n)T n (x) . (2.17) 

n=l 

Clearly (recalling (|2.6p ). p c (#) = (<Siv/c) for all G [— 7r,7r] . Since Tj^jv is exact for all 
polynomials of degree up to N, we have, using Theorem 12.11 



|4(/)-4,jv(/)| = l4(/-p)-4,iv(/-p)| 

= ^-J (^J \\(i-SN)f c \\ Hm . 

Then, using Lemma 12.31 

/ 1 \ r / i \ w»-p(r) 

|Jfc(/)-/fc,Jv(/)| < C^atM (^J \\(f (m) )c\\ H o (2.18) 

and the result follows from (|2.15p . □ 
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2.2 Integrals over [a, b] 

Now we consider the integral (|2,ip for general [a, b]. To apply the FCC quadrature, we first 
transform the integral using the following linear change of variables: 



x = c + ht, t€ [-1, 1], 
Then we may write 



b + a , , b - a 

where c := , and h := 

2 ' 2 



I^(f) = hexp(ikc)I~ k (f), 
where k = hk and / is the function on [—1,1]: 

f(t) = f( c + ht) , te [-1,1] . 



(2.19) 
(2.20) 

(2.21) 



Then we apply the quadrature rule (|2.2p - (|2,4p to the integral on the right-hand side of (|2.20|) 
to obtain the approximation 



I l «§(f):=he X p(ikc)I kN (f) » if •"'(/)• 
The following theorem is the corresponding extension of Theorem | 

Theorem 2.6 Let r G [0, 2] and < m < N + 1. Then, when \ f\H™[a,b] < °°> we have 



(2.22) 



t b \f)-ltN(f) 



a,b] 



< at 



1 A; J hm+1 ~ r ( ]y ) l/lfl&M] 



(2.23) 



Proof. From (I2.20|) . (I2.22D and then Theorem 12.51 we obtain 



/^ ] (/)-«(/) = h I~ k (f~) - I kN (f~) 
< a' mN h(V\ | — 



fc7 Viv 



m— p(r) 



1 ^ 



l/~k 

-p(r) 



(2.24) 



Now = h m /M (c + /it), and 



so 



and the result follows. 



2m 



|/( m ) (c + /it)| 



dt = h? m \f\l s[aA 



(2.25) 
□ 



The most important use of this theorem will be for the case when N is fixed and convergence 
is obtained by letting h — >■ (as arises when composite versions of the FCC rule are used). For 
this case we have the following corollary, which is obtained using Theorem 12.61 with m = N + 1. 

Corollary 2.7 Letr G [0,2]. For each N > 1, there exists a constant c N = Ca' N+1N (l/N) N+1 - p ( r \ 
such that 



< CAT 



(1.Y h N+2 ~ r max 

\kj x€[a,b] 



(2.26) 



when f G C N+1 [a,b}. 
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3 Composite Clenshaw-Curtis Rules 



In this section, we will consider the computation of (/), where / is allowed to have an 
algebraic or logarithmic singularity in [a, b). To control the length of the paper, we restrict to 
functions / which are not continuously differentiable. Singularities in higher derivatives can be 
treated in an analogous way. 

Without loss of generality we set [a,b] = [0, 1] and assume that the only singularity occurs 
at the origin. The case of a finite number of singularities on [a, b] can be treated by splitting 
[a, b] up into subintervals, each with only one singularity at an end point, and them mapping 
each interval onto [0, 1] in an obvious affine way. Hence, for /? E (0, 1) and m > 1, we introduce 



\ v \ |m,/3 := m ax < sup sup 
U-e[o,i] xe(o,i] 



j = l,...,m}. (3.1) 



We denote by C™[0, 1] the space of all functions v E C[0, 1] such that ||u|| m ,/8 < °°- Similarly, 
for f3 E (—1,0) we define 



\v\ \m,p '■= m a x \ sup 

I x6(0,l] 



x V-P) v V)(x) , j = 0,...,m} , (3.2) 



and choose C™[0, 1] to be the space of all v E C(0, 1] such that |H| m ,/3 < oo. Finally, we cover 
the case of logarithmic singularities via the norm 

IMIm.,o := max { sup |(| logx\ + l) _1 v(3;)|, sup x^v^\x) , j = l,...,m> (3.3) 

l_a;6[0,l] x6[0,l] J 

and introduce the associated space C^[0, 1]. Note that C^[0, 1] C L^O, 1] for all f3 E (—1,1). 
3.1 The composite algorithm 

When / E C™[0, 1], for f3 E (—1,1), our strategy for computing is to apply the FCC 

rule in a composite fashion on a mesh graded towards the singularity. With the right choice 
of mesh the error of the quadrature can then be made to satisfy a uniform error estimate on 
subintervals and be small overall. Let us recall the classical graded mesh 

U M , q := Ixj := (j^\ q : j = 0, 1, . . . , m\ , (3.4) 

where q > 1 is the grading parameter to be chosen. This mesh - originally proposed in |19j - is 
well-known to give optimal approximation of functions with singularities by fixed order piecewise 
polynomials. An application to quadrature was given in [7]. This paper contains an extension 
of these results to the computation of oscillatory integrals with singularities. Writing 

M 



I f\ f ) = i^\f) + Y^it j - UXi \f), 

i=2 

we approximate each term in the sum on the right-hand side by applying the FCC rule as 
defined in (|2.22j) . The strategy for approximating the first term on the right-hand side depends 
on whether /3 < or (3 > 0. Precisely we define the approximation 

h (/) -"io, if pe (-1,0]. (3 ' 5) 



9 



Note that for /3 G (0, 1), 

i£r i] (f)= 



I 1 (q[xoMj\ ( x ) ew (±kx)dx, if xik > 1 
Q [ *°' Xl] (/exp(ifc-))) (x)dx, if x x k < 1 



fti 



■>-■<) 



where Q^' Xl 'f is the linear function interpolating (xq, f(xo)) and (x\, f(x\)). (To obtain this 
formula, recall that from §Z7M . if 3 ^ (/) = he-xp(ikc)I- k x (f) where A; = fcci/2, and recall (pT2|) 
and (|2.4p ). The composite quadrature rule is 



A I 



1 k,N,M,q\J I ■— 1 k U I "T Ui' 

J'=2 

The corresponding error may then be bounded by 
where 



(3.6) 



i=2 



ei 



/f°' :El] (/)-^ ' a:il (/) ) 



and 



^ 1,Xjl (/)-^ 1 ' Xj] (/), for J=2,...,M. 



(3.7) 



(3.8) 
(3.9) 



In the following two sections, we derive results which will help us estimate \e~i\. These are 
subsequently used to estimate the total error f fc jv,Af,g(/) i n Theorem 13.61 

3.2 Estimates on the size of the integrals 

In the following two lemmas, we analyse the integrals I k °' £ \f), where / G Ca[0, 1], making 
explicit the rate of decay as both e — > and k — > oo. 

Lemma 3.1 For any (3 G (—1,0) and any f G Cg[0, 1] f/iere exists Cp > such that fore G (0,1] 



i4 0,£l (/)i < c^-Q 



l,/3 



(3.10) 



where s G [0, 1+/3]. Furthermore, for any f G Cg[0, 1] i/iere exzsis Co > suc/i that for e G (0, 1], 
we /iai>e 



/f eJ (/)l < G)(e + £|loge|) J 



i_ s / 1 + log k 
k 



1,0 



(3.11) 



where s G [0, 1] . 



Proof. We will prove the lemma for the case when /3 G (—1, 0). The case when /3 = follows 
similarly. First note that for all e G (0, 1], we have 



i4°' £l (/)i < 

We show now 



x^ dx 



1 



0.0 



1 + /3 



£ 1+f3 \\f\\o,p < 



1 



i4 0,£l (/)i < 



2 + 



(3.12) 



(3.13) 
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and the result follows by interpolation of (|3.12p and (|3.13p . 

To obtain (|3.13p . note first that it follows trivially from ()3. 12|) if ek < 1. Therefore, let us 
assume that e > l/k, and write 



#.-J( /) = J f.V*J (/) + i iVM (/) . 



Now again from ()3. 12j) 



r [0,l/fc] 



(/)l< 



1 (\ 



On the other hand, integration by parts yields 



ll ,k ' £ \f) = \\f{x)e W (ikx)V ] -i 

IK L J x=l/k lk 



0,P- 



f'(x) exp(ifcx) dx. 



l/k 



Thus 



l/f/M^i < 



+ |/(1/*)| + / |/'(x)|dx 
l/fc 



Now for any a; > 0, < anl/Hi,/? and also 



l/k 



\f'(x)\dx < / x^&xWfW^ < 



l/k 



Thus 



1 + ¥\ 



and, since < (l/k) 13 we obtain 

|/f /M (/)l < 2 N 



101 

Substituting (l3~T6|) and (l3J5|) into {gJSj) we obtain 

+ 2(1-|/3| 2 ) 



(1 



1 \ 1+/3 



(3.14) 
(3.15) 



(3.16) 



thus proving (|3,13p . 



□ 



Lemma 3.2 For any (3 G (0,1) and for any f G Cl[0, 1] there exists Cp > such that for 
e G (0, 1] and s G [0, 1] we have 



\if e \f)\<c^{^j 



(3.17) 



Moreover, if s G [1, 1 + /?], 



i4°' £J (/)i < 



1/(0)1 + 



+ C p e 



l+P-8 



1,8-1 



2,/3- 



(3.18) 
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Proof. First note that for e G (0, 1], 



i? ,Bl (/)|< ell/Ik*. 



On the other hand, integration by parts yields 



t £ \f) 



< 



k 



1/(0)1 + 1/(6)1 



+ 



f'(x) exp(ifcx) dx 



f'(x) exp(i/cx) da; 



Also, we have 



f'(x) exp(ikx) dx 







< 


[ x^dx 




Jo 



(3.19) 

(3.20) 
(3.21) 

(3.22) 



and substitution of (I3.22|) into (|3.2ip yields 



Interpolation of this with (|3.19p yields (|3.17p . 

To obtain (|3.18p . we also note /' G Ci-JO, 1] so by Lemma I3TT1 we also have 



f'(x) exp(ifcr) dx 



2,0 



for all s' G [0, f}\. Substituting this into (|3.20p and putting s = 1 + s', we obtain the result. □ 

In the next lemma we shall verify the sharpness of the estimates in Lemmas 13. HI3.2l as k — > oo. 
This result concerns the family of functions 



Mi 



/?G (-1,0) U (0,1), 



log x , fi = 0. 



(3.23) 



Lemma 3.3 For all (3 G (—1, 1), there exists a constant Ap > such that, for all k sufficiently 
large. 



fc min{imi}| 7 [0,l] (//3) | > Afjf when 

l4°' 1] (/o)l > U 



log k 



(3.24) 
(3.25) 



Proof Let us consider first (3 G (-1, 0) U (0, 1). Using (3.761.1), (3.761. 6)], we obtain 

1 



1 + 



iFi(l + 0,2 + 0,ifc), 



(3.26) 



where \F\(a, b, z) denotes the confluent hypergeometric function (also called Kummer's function 
and denoted M(a,b,z) in [lj Section 13]). At large values of \z\, with — it/2 < argz < 37r/2, for 
fixed a and b, the function ±Fi(a,b, z) has the following asymptotics, see [H (13.5.1)]: 



iFx(a,b,z) = T(b) 



T(b-a) 



+ 



e z z a-b 



(i + oCM- 1 )). 
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Then, from (|3.26[) and since T(l + z) = zT(z), we obtain, 
Therefore, for k sufficiently large, we have, for /3 G (—1,0), 

and for /3 G (0,1), 

\t 1] m > Kk-i-ni+w-p) > \k-\ 

These prove the estimate (|3.24|) . 

To verify (|335jh we use the formulae [U (4.381.1), (4.381.2)] to obtain 

4°'%o) = ~\ (| + i7 + itogfc) + £(<*(*) + isi(fe)) 

where si and ci are the sine and cosine integral functions: 

... 7r /" x sin i ... rcost-l , 

si(x) := — — + / — — dt, ci(x) := 7 + logx + / dt, 

* Jo t Jo t 

and 7 ~ 0.5772 is the Euler-Macheroti constant. (Note that, in the notation of pQ, si(x) = 
— 7r/2 + Si(x) and ci(x) = Ci(x).) Then using the asymptotics for large arguments of Si and Ci 
in [U (5.2.34), (5.2.35)] and pQ (5.2.8), (5.2.9)], we deduce that 



lim si(fc) =0(l/k), lim ci(fc) = 0(l/k). 

k— >00 fc— >QO 



Thus 



^(/.) = -i(f + i 7 + n„ g(: ) +0 (J,). 

Thus, for all k sufficiently large 

1 k ~ ~2k~' 

proving the estimate (|3.25[) . □ 



3.3 The total error for the composite Filon-Clenshaw- Curtis method 

In Theorem 13.61 below we use (|3.7p to estimate the total error of the composite FCC rule. The 
first contribution \e\\ is estimated either by a direct application of Lemma l3.1l (when (3 G (—1,0]), 
or via an integration by parts argument (when /3 G (0, 1)). This is done in Lemma 13.51 below, 
but first the remaining sum on the right-hand side is estimated in the following lemma. Since 
the proof uses fairly classical graded mesh arguments we shall be brief. 

Lemma 3.4 Let f G C^ +1 [0, 1], /3 G (-1,1), letr>0 and choose 

q > (JV + 1 - r)/(/3 + 1 - r), for r < 1 + /3. (3.27) 

Then there exists a constant C which depends on N , /3 and q such that 

£N * c {k) (m) H/IIw ( 3 - 28 ) 
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Proof. In the proof we let C denote a generic constant which may depend on N, (3 and q. 
By Corollary 12. 7| and denoting hj = (xj — Xj-i)/2, we have 



M , x r M 

Y\ e \ < C ~) Yh? +2 ~ r max 



(x) 



■ t \ r I -M 

i=2 



JV+2-r 



iV+1,/3 



A simple application of the mean-value theorem shows that 



1 fj-1 



h, < C- 
3 ~ M M 



q-l 



, for j = 2,...,M, 



and hence 



o J- 1 - \M) 

where a = q((3 + 1 — r) — (N + 2 — r), and so 



-, \ N+2-r / ■ -. \ a 



M 



M / 1 \ r / -I \ JV+l-r | M 



J'=2 



^ M V M 



i=2 



JV+1,/3 ' 



(3.29) 



(3.30) 



The result follows since the factor in braces in (|3,30p is a Riemann sum for the integral x a dx 
- which is finite, since the hypothesis of the lemma ensures that a > — 1. □ 



Lemma 3.5 Under the same hypothesis as Lemma \3.4\ there exists a constant C which depends 
on (3 and q such that, for N > 1 

/l\ r / 1 \ N+1 ~ r { ' when P G (- 1 ' ) u (°: !); 

\k) \Mj \ (l + logfc) r (logM) 1 - ||/|| 20 , when = 0. 

Proof. Throughout we use the fact that xq = 0. Consider first j3 G (0, 1) and note that 



\(Q?' xl] f)'\ 



f(xi) - Ho) 



l 

< — 

Xl 



Xl 



1 f X1 I 

— / f'(x)dx 
xi Jo 



x' dx 



1 f: 



i r xi 

< _ / \f'( x )\dx 
xi Jo 



l,/3 • 



Moreover, for any t G [0, xi] and any / G C«[0, 1], we have 



| /W _ Q f^] /(t) | 



{f-Qf Xl] f)'{x)\dx 

10 

< n\f{x)\dx+ xl \(Q^fy\< -Awn^ 



JO,Xl] ,\/| 











Then, from ([33]), we have from (|332|) . when / G C$[0, 1], 

^ (/(x)-Ql ' Xll /(^))exp(ifcx)dx < |a> 



1„8 



(3.31) 



(3.32) 



(3.33) 
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On the other hand, integrating the formula for ej by parts, we obtain 



ei 



1 f Xl 
ik 



f(x)-(Q [ ?' Xl] f)(x)) exp(iA;x)dx . 



(3.34) 



Since /' G Cfl_i; to treat the first term on the right-hand side of (13.34p we can use Lemma 13.1 
with [3 replaced by [3 — 1 and s chosen to be /3, thus obtaining 



f'(x) exp(i/cx)dx 



< Ca - 



2,13 



(3.35) 



Moreover, to treat the second term in (|3.34p . since (Q^i' Xli f)' is constant, we have by (|3.3ip 



{Q [ °' Xl] f)' exp(ikx)dx 



exp(iA;x)dx 



< 



Hence combining f)3.35|) and (|3.36p with (|3.34[) . we obtain 



\ei\ < C, 



+ 



1 ~ C>13 I k 



P 



1 + (a?ifc) 



/3Jfe 



l,/3 



(3.36) 



A;/ A; 2 

Hence if xife > 1 we can interpolate (|3.37p and (|3.33j) . to deduce that 



(3.37) 



JV+l- 



2,0 



(3.38) 



for any r G [0, 1 + /?]. 

On the other hand, if kx\ < 1, and defining := /(x) exp(i/cx), (I3.32p yields 



e i 



2 / 1 



{fk(x)-Qf Xl] f k {x))dx 



< -= t- a; 



„l+/3-r 
'1 



b/3 



2/iy/p * 



(3.39) 



and (|Ofl) . (^391) prove the result for /3 G (0, 1) . 

For /3 G (— 1, 0] note that the integral over the first subinterval is approximated by zero, and 
so the result follows readily from Lemma 1331 (with e = x\ and s = r). □ 



The proof of the following result now follows directly from Lemmas 13.41 and 



Theorem 3.6 Under the same hypothesis as Lemma \3.4\ there exists a constant C which de- 
pends on N , (3 and q such that 



Ek,N,M,q(f) < C 



1 

17 



. when /? G (-1,0) U (0,1), 

(1 -f logfc) r (log Mf~ r \\f\\ N+lfl , when = 0. 
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4 Nonlinear oscillators 



In this section we return to the integral l[ a '^(/, g) defined in (jl.ip . and consider a general 
nonlinear g. We will assume for simplicity that g E C°°[a,b]. For less smooth g the arguments 
will be analogous but the exposition would be more technical. Our methods will be based on 
the change of variable 

T = g(x). (4.1) 

If g has no stationary points (i.e. g' does not vanish), then g" 1 € C°°[a, b], {g~ 1 )' = l/(g' ° 9 l ) 
and 

lt b] (f,9) = A 9(a) ' 9ib)] (F), With F = (f o g-l)\(g-ly\ = (f o g-l)\(g> o g-l)\-\ (4.2) 

Now, assuming also that / € C°°[a, b], then (|4.2p can be computed using the FCC rules described 
in $21 with the additional cost being the evaluation of the inverse function g~ l at the quadrature 
points. Moreover if / has singularities then these induce singularities in F and the composite 
FCC rules in $3] could be used instead. (Here we assume implicitly that g(a) > g(b). If g{b) > 
g{a) then the integral can be transferred to one over the interval [g(b),g(a)\ by a simple affine 
change of variables. We make similar implicit assumptions below.) 

If now g has a stationary point at one or more £ E [a, b], the transformation (|4.ip may still 
be applied, but singularities appear in F at the points <?(£)• To describe these, we may, without 
loss of generality, consider a single stationary point £ E [a, b] of order n > 1 with property 

g'(0=g"(0 = --- = g {n \0 = 0, 9 (n+1) (0>0 and g'{x) + 0, for x E [a,b] \ {£}. (4.3) 

Then g is monotone on each of the intervals [a, £) and (£, 6]. The change of variables (14. ip can 
be applied on each interval separately to obtain 

lt M Uig) = [I +/ )F(r)exp(ifcr)dr (4.4) 

with F as in (|4.2[) . (One of the integrals in (|4.4p is interpreted as void if £ = a or 6.) The regu- 
larity of the resulting function F is summarised in the following theorem. Here and throughout 
the rest of this section we use the convenient notation a = l/(n + 1). 

Theorem 4.1 Assume that f and g are in C°°[a,b] and that g has a single stationary point £ 
of order n as in (|4.3p . Then for each p £ No := N U {0}, there exists C p > such that 

F<p\t) < C p \r - 5 (0r- p - 1 , TG[g( a ), 5 (0)U (5(0,5(6)] • (4-5) 

The proof of Theorem 14.11 requires Lemma [4.2l (below) and both these results require the Faa 
di Bruno formula for the derivatives of the composition of two univariate functions. For p G No, 
let denote the pth derivative of any sufficiently differentiable function /. Then, if <f>, tjj are 
suitably smooth functions and the composition (p o ip is well-defined, we have the formula 

(^) (p) =E^(> (H) °^ (n w i - ( 4 - e ) 




where the sum in (|4.6p is over all multiindices m G (Nd) p which satisfy mi + 2m-2 + . . . +pm p = p. 
Moreover \m\ = mi + m,2 + . . . + m p and ml = m^.m^}. . . . m p \. A suitable reference is ( [201 
Theorem 2]). 
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Lemma 4.2 Assume g G C°°[a,b] and that g has a single stationary point £ of order n as in 
T3|) , Then, for all p£N, there exists a constant C p > swc/i £/ia£ 

GT 1 ) 60 ^)! < C p \r-g(0\ a ~ p . (4.7) 

Proof. Without loss of generality we assume £ < 6, that <? is increasing in [£, 6] and we consider 
the case r G (£, 6] only. (The case r G [a, £) is completely analogous.) By Taylor's theorem with 
integral remainder and (|4.3p . 



5 (x) = ff (£) + (x- 09(0 + ■■■ + ^P^9 (n) (0 + R*(x) = g(0 + Rz(x), (4.8) 
for all x G (£,6], where 

Rt(x) = — I* {x-t) n g( n+l \t)dt. (4.9) 



With the change of variables t \— > y = (t — £)/(x ~ i R <|4.9|) . we obtain 
i? 5 (x) = (x-£r +1 T 5 (x), where T f (x) = i | (1 - y) V n+1) (£ + - 0) dy. (4.10) 

Then, for all x G (f , 6], # 5 (x) > and T^x) > 0. Also, since g G C°°[a, 6], we have T ? G C°°[£, b] 
andT € (e) = ^ 5 (™+ 1 )(0>0. 

Now, recall q = l/(n + 1) and define 

he(z) = (g(x)-g(0) a = (Rdx)T = (mx)) a (x-0- (4.11) 

Then ^ G C°°(£, b] and, for all x G (£,6] , /i^(x) = a(g(x) - giClT" 1 g' '(x)> . Moreover since 
also h'g(£) = (T^(£)) Q > it follows that h'^ is positive valued on [£,6] and so : [£,6] — > R is 
invertible and (/if) -1 G C°°[^(£), ^(6)]. Thus, inserting (03) (and x = ff _1 (r)) into (|4TTT|) . we 

^(t) = x = ^((t-^))-). (4.12) 

To prove the estimates (|4,7p we now apply the Faa di Bruno formula (|4,6p with (f> = h^ 1 and 
V'( r ) = ( r — 9{0) a to obtain derivatives of (|4.12j) . Consider any term in the resulting sum (14.61) . 
Since /i^" 1 is smooth, the first factor in round brackets is bounded, while the second factor in 
round brackets can be estimated by 

| T _ ^^|(a-l)mi+(a-2)m2+...+(Qt-p)m P) (4.13) 

times a constant. Recalling the remarks following (14. 6j) . the index in (|4.13j) is a|m| — p > a — p, 
so (02)) follows. □ 



Proof of Theorem \4-l\ Again, without loss of generality, we work with r G (#(£), g(b)]. We first 
observe that since g' is one-signed, so is (g" 1 )'. Thus we can write F = ±(/ o g~ l ) (g^ 1 ) and 
hence, by the Leibnitz rule, F^ is a linear combination of terms of the form 

(f o g^f (g- l ) {p - l+l) , 1 = 0,..., p. (4.14) 

Referring again to formula (|4.6p . and recalling that / is smooth, the first term in (|4.14p may be 
estimated by a constant times 

I 

i=i 
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where mi + Im^ + • • • Imi = I. Using the same argument as in the proof of Lemma 14.21 this 
product has the estimate |r — g{£)\ a ~ l (modulo a constant factor). 

Now, returning to the products (|4.14p . we see that for I ^ each of these can be estimated 
(modulo a constant factor) by 

|r - g(Or l \r ~ g(Or P+l - 1 = \r~ <7(0| 2a " P-1 ■ 

However, when I = the bound is |r — g(£,)\ a ~ p ~ l an d since a > 0, the result (|4.5p follows. 
□ 

4.1 Accurate implementation 

Now let us return to the computation of (|4.2p . Under the assumption of Theorem 14. 1\ we write 

I [ k' ] (f,9) = / +/ F(r)exp(iA;r)dr . (4.15) 

\Jg(a) Jg(0 J 

Each of these two integrals can be transformed in an afhne way to an integral over [0, 1] so that 
the singularity is placed at the origin. The composite FCC algorithm given in £ j3.1l can then be 
applied, with error estimates given by Theorem 13.61 For example, consider the second integral 
in (|4.15p . Under the change of variable r = <?(£) + cx where c = g(b) — <?(£) and x S [0, 1] this 
becomes 

4° C ' 1] (F) where F(x) = cexp(i%(£))F(s<£) + cx) , (4.16) 

and, by TheoremHIEJ F E Cp[0, 1], with p = a -l = -n/(n + l) G (-1,0). 

In the implementation of the composite FCC rules for (|4.16p some care must be taken to 
accurately evaluate the integrand F(g(£) + cx) at very small arguments x (as arise in the case of 
finely graded meshes). This is a delicate matter since if g{£) 3> cx, rounding error may pollute 
the direct calculation of <?(£) + cx, in turn making F(x) inaccurate. To solve this problem, recall 
that F is defined in (|4.2p in terms of the composition of smooth functions / and g' with g^ 1 . Our 
task is therefore reduced to devising an accurate evaluation of the quantity x := g~ 1 (g(C) + e ) 
for small e. 

The required x is then a solution to the equation g(x) — = e and, recalling the proof of 
Lemma H21 we see that this is in turn equivalent to (T^(x)) a (x — £) = e a . Thus x solves the 
nonlinear parameter dependent problem 

G(x,e) := m(x)) a (x-0-e a =0 . 

Since T^(£) > (see also the proof of Lemma [472]) . we have G(£,0) = / G x (£, 0) and so the 
Implicit Function Theorem implies that, near e = 0, x is a smooth function of e a and there 
exists a constant C\ so that \x — £| < C\e a , for small enough e. Moreover x is also a solution to 
the fixed point problem 

*= ?+ Gik)" =: " (a;) - 

Since (£) > and is smooth in a neighbourhood of £, it is easy to see that H is Lipschitz in 
a ball centred on £ and its Lipschitz constant is C^e" for some constant C2. So for small enough 
e, x is the unique fixed point of H and fixed point iteration converges. This suggests that for e 
small, a suitable approximation to x can be chosen as x := H(£), with error 

\x-x\ = \H(0-H(x)\ <C 2 e a \i-x\ = 0(e 2a ) . 

The approximation x to 5'~ 1 (5'(C) + e) when e is small is used in the computations in §5.21 
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5 Numerical Experiments 

In this section, we first carry out some numerical experiments which illustrate the convergence 
estimates of Theorem 13.61 using computations of the model integral with linear oscillator: 



for various /5 € (—1,1), with fp defined in (|3.23p . Then we compute a model problem with a 
nonlinear oscillator motivated by the implementation of hybrid numerical-asymptotic boundary 
integral methods in high-frequency scattering. 

5.1 Linear oscillator 
Experiment 1 

Our first set of experiments studies the case M — > oo, for fixed k. From Theorem 13.61 with 
r = 0, we see that in this case the composite FCC rule for (|5,ip should converge with order 
0(M"( Ar+1 )) as M — ^ oo provided q > (N+l)/(/3+l) for (3 / 0. When j3 = an additional factor 
of logM appears in the estimate. To illustrate this result we compute the errors -E/c.^m^C/^) 
for k = 1000 with various N and q = (N + l)/(/3 + 1) + 0.1 as M increases. The exact value 
of (|5.ip can be computed analytically and so the errors can be found exactly. The results for 
the three values f3 = 1/2,0,-1/4 are given in the three sub-tables in Table [TJ The columns 
headed "error" contain the values of ^A;,iV,M,g(/^) while the columns headed "ratio" contain 
the empirical convergence rates with respect to M computed by extrapolation The expected 
convergence rate is N + 1 (modulo a log factor when f3 = 0) and this is given in the row marked 
"expected ratio". In all cases the empirical convergence rate is close to the predicted rate, 
except when the error has almost reached machine precision in which case, naturally, rather 
unsteady empirical convergence rates are obtained. It is worth noting that in this computation 
some of the subintervals in the composite rule are very small, in fact with N = 8 and M = 64 
the smallest subinterval of the mesh is of size about 10 -34 . Nevertheless the algorithm appears 
to show no instability and converges to machine precision as M increases. 

Experiment 2 

Here we fix M = 10, N = 3 and q = 12 and we study convergence as k increases, for various j3. 
In Table O the columns headed "ratio" contain the empirical convergence rates with respect 
to k computed by extrapolation. From Theorem 13.61 we see that the composite FCC rule for 
(|5.ip should converge with order 0{k~ r ) as k — > oo where r < (q(f3 + 1) — N — l)/(q — 1). In 
the row marked "best expected ratio" this upper bound on r is given for each f3, using our 
choice of N,q. We see from Table [2] that the empirical convergence rate as k increases for f3 > 
is close to the theoretically predicted best rate of convergence. When (3 < 0, the empirical rates 
are a bit slower that the theoretical best rate. 

Experiment 3 

In Table [3] we study the computation of (|5.1|) with //j(x) = logic, for M = 12 and N = 3 as k 
increases for various q. For each value of q and N, the error of the composite FCC rule should 
converge with order 0(k~ r ), with r < (q — N — l)/(q — 1). The row marked "best expected 
ratio" contains the upper bound on r while the columns marked "ratio" contain the empirical 
convergence rates . The table shows that when when q = 12 the empirical convergence rate is 




(5.1) 
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= 1/2 
expected ratio 


TV = 4 
5 


TV = 6 
7 


TV = 8 
9 


M 


error 


ratio 


error 


ratio 


error 


ratio 


8 
16 
32 
64 


4.3e-006 
9.5e-008 
2.9e-009 
8.1e-011 


5.49 
5.03 
5.17 


5.2e-008 
5.7e-010 
2.0e-012 
2.3e-014 


6.50 
8.13 
6.50 


1.7e-009 
6.6e-012 
1.0e-014 
1.3e-016 


8.06 
9.30 
6.37 




/3 = 
expected ratio 


TV = 4 
5 
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4.10 
7.15 
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7.62 
6.43 
7.88 


6.0e-006 
2.0e-008 
l.le-011 
2.9e-014 


8.22 
10.91 
8.51 



Table 1: Numerical Results for Experiment 1 



close the the theoretical best rate. When q = 4, the empirical convergence rate is even better 
than the best expected rate which in this case indicates that no convergence should be observed 
at all. On the other hand, when q = 8 and q = 16 rather unsteady convergence rates are 
obtained. However, when q = 8 as k increases the empirical rates become bounded by the best 
expected rate, while for q = 16 the empirical rates are either bounded by or are slightly better 
than the the best expected rate. 

Experiment 4 

Here we illustrate the power of the composite FCC rule compared to the non-composite version 
for computing f)5. 1 1) when (3 = 1/2. The parameter TV for the non-composite FCC rule takes 
values Ni = {24 x 2 l ,i = 0, 1,2,3}, while for the composite rule, we fix parameters q = 12 and 
M = 6 and take iVj = {4 x 2*, i = 0, 1, 2, 3}. The total number of function evaluations in both 
cases is therefore the same. The superiority of the composite version is clearly seen. 

5.2 An example from boundary integral methods in high-frequency scatter- 
ing 

Finally, we describe an application to the computation of acoustic scattering at high frequency 
by numerical-asymptotic methods. When an incident plane wave exp(iA;a;.d) is scattered by a 
smooth convex sound-soft obstacle with boundary T, the scattered field can be computed by 
solving the integral equation 

J ^\k\K - y\)v(y) ds(y) = exp(i£:x.d) , x G T. 

where ^q 1 ^ is the Hankel function of the first kind of order 0. (In fact this problem is not 
well-posed for all values of k and in practice a related "combined potential" formulation is 
used. However this formulation illustrates the essential quadrature challenge which arises in all 
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Table 2: Numerical Results for Experiment 2 
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Table 3: Numerical results for Experiment 3 
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Table 4: Numerical results for Experiment 4 
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formulations. ) Even for moderate values of k is useful to apply the "physical optics approxi- 
mation" which amounts to writing v(y) = V(y) exp(ifey.d). and computing the less oscillatory 
component V rather than the highly oscillatory v - see [3], [2]. 

Using the fact that H^ikr) exp(— ikr) is a non-oscillatory function, smooth for r > but with 
a logarithmic singularity at r = and introducing a smooth parameterization x : [0, 2tt] — > T, 
the problem above can be reformulated, for s £ [0, 2ir] as 



[ 2W M k (s,t)exp(ik^ [s] (t))V(t)dt= 1, * [a] (t) := |x(s) - x(i)| - (x(a) - x(t)).d . (5.2) 
J o 

The function Mf~(s,t) is non-oscillating and smooth except as i = s where a logarithmic singu- 
larity occurs. (More details are in [T3] .) It can be proved that when s is chosen so that x(s) 
is in the "illuminated" part of T (where the incident waves hits the obstacle), there is only one 
stationary point, i.e., there exists a unique t so that (^[ s ]) (t) = 0. Moreover, x(t) is a point in 
the shadow part. Conversely, if x(s) lies in the shadow, we have three stationary points, with 
two of them in the shadow and one in the illuminated part. Thus (|5.2|) is a good example for 
application of the methods of $31 

As an illustration we compute the integral in (15. 2j) when T is the unit circle, x(s) = 
(cos s, sins) and set d = (1,0) so that 



*w(t) 



sin 



s 



cos s + COS t. 



The integral is computed using the following strategy. First, [0, 2ir] is divided into subintervals so 
that each of them contains at most one point which is either the singular point s or a stationary 
point. Next, any subinterval of length greater than one is split into two subintervals of equal 
length and this process is continued until we obtain a subdivision of [0, 2tt] in say J subintervals 
of length smaller than 1. This is done to avoid working with graded meshes on relatively long 
intervals. Then, with the change of variable r = \I/r s i(i), we have to approximate the integral 



/ W 3 F s (r)exp(ifcr)dT, F s (r) := M k (s,^ (t))V^,Ut)) } j = l,...,J. 



Given L a positive integer, we compute the approximation of these integrals as follows: If [aj, bj] 
does not contain s or a stationary point, the function F s is smooth, and therefore we can apply 
the simple Filon Clenshaw-Curtis rule with min{L + 1, 129} points. If s S [aj, bj], after an affine 
change of variables, F s belongs to Cq°[0, 1]. On the other hand, if either aj or bj is a stationary 
point of order n > 1, F eCf^// +] \ [0, 1] cf. Theorem 14.11 In both cases, we introduce meshes 
appropriately graded towards the singularity, according Theorem 13.61 with L subintervals and 
apply the composite FCC rule with iV + 1 points. In our experiment we have taken V = 1 
and s = Sir/A, which corresponds to a point in the illuminated part. The only stationary point 
t = 23vr/24 is of order 1. 

In Table [5] we display the results obtained with N = 6, q = (N + 1)/(1 + (3) + 1/10, which 
corresponds to r = in Theorem 13.61 Here (3 = or —1/2 depending on the singularity of the 
integrand. From Theorem 13.61 we can expect that the error decreases with as 0(L~ 7 ) (but there 
is no predicted decay with respect to k in this case) . In Table [6l we give results for N = 4 
and q = (N + 1)/(1 + j3 - 1/4) + 1/10, leading to the error estimate 0(k~ 1/A L~ 4 - 75 ). Although 
the convergence in Table [6] is irregular, in contrast to Table [5j the error decreases with k as 
predicted by the theory. 

In this experiment over 90% of the CPU time is spent in computing the change of variable, 
since any evaluation in r requires the solution of a non-linear equation. This is carried out in our 
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implementation using the f zero of Matlab. In the case that (|5.2p has to be computed for many 
different functions V, for instance in assembling the matrix of a Boundary Element Method, 
these calculations have to be done only once, resulting in a speed up of the the method. We run 
our programs in a modest laptop, a Core 2 Duo with 4Gb of Ram memory, and show in Table 
[7] the CPU time required to construct some columns of Table [5l We clearly see that the cost is 
linear in L. 
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Table 5: Results for N = 6, M = L and q = {N + 1) for intervals with log singularities and 
q = (N + !)/(! — 1/2) for integrals having a stationary point (so r = 0). 
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Table 6: Results for N = 4, M = L and q = (N + 1)/(1 — 1/4) for intervals with log singularities 
and q = (N + 1)/(1 — 1/2 — 1/4) for intervals having a stationary point (so r = 1/4). 
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Table 7: CPU time consumed for computing Tabled] 
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